#include "../DWMRICore/Fiber.h"
#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/DTI.h"

#include "CUDA_ComputeStructureTensor.h"

#include <omp.h>

char *fiber_pathname = NULL;
char *dst_pathname = NULL;
char *tempPathname = NULL;
char *fa_pathname = NULL;

char **confPathname = NULL;
char **vertPathname = NULL;
int w, h, d, scale;
float stepSize, maxLength;


CDTI *pStructureTensor = NULL;


float fiber_length = 0.0f;
int kernel_size = 0;
int mode = 0;		/* 0 for dti, 1 for hardi */

enum ST_METHOD {
	METHOD_NONE,
	METHOD_POSITION,
	METHOD_FA,
	METHOD_EV
};

ST_METHOD method = METHOD_NONE;

extern void ComputeStructureTensorMethodPosition(const float fiber_length, const float maxLength, 
	const float stepSize, const int w, const int h, const int d, const int scale, const int mode, 
	char **confPathname, char **vertPathname, char *tempPathname);

extern void ComputeStructureTensorMethodFA(const float fiber_length, const float maxLength, 
	const float stepSize, const int w, const int h, const int d, const int scale, const int mode, 
	char **confPathname, char **vertPathname, char *tempPathname, char *faPathname);

extern void ComputeStructureTensorMethodEV(const float fiber_length, const float maxLength, 
	const float stepSize, const int w, const int h, const int d, const int scale, const int mode, 
	char **confPathname, char **vertPathname, char *tempPathname);


void PrintUsage()
{
	printf("Usage: ComputeStructureTensor -fiber <configure_file> -dst <dst_file> -length <fiber_length> -kernel <kernel_size> -mode <mode> -method <method> [-temp <temp_file> -fa <fa_file> -md <md_file>]\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-fiber") == 0) {
			fiber_pathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			dst_pathname = argv[++i];
		} else if (strcmp(argv[i], "-fa") == 0) {
			fa_pathname = argv[++i];
		} else if (strcmp(argv[i], "-length") == 0) {
			fiber_length = atof(argv[++i]);
		} else if (strcmp(argv[i], "-kernel") == 0) {
			kernel_size = atoi(argv[++i]);
		} else if (strcmp(argv[i], "-temp") == 0) {
			tempPathname = argv[++i];
		} else if (strcmp(argv[i], "-mode") == 0) {
			i++;
			if (argv[i][0] == 'D' || argv[i][0] == 'd') {
				mode = 0;
			} else if (argv[i][0] == 'H' || argv[i][0] == 'h') {
				mode = 1;
			} else {
				return -1;
			}
			
		} else if (strcmp(argv[i], "-method") == 0) {
			i++;
			if (argv[i][0] == 'P' || argv[i][0] == 'p') {
				method = METHOD_POSITION;
			} else if (argv[i][0] == 'F' || argv[i][0] == 'f') {
				method = METHOD_FA;
			} else if (argv[i][0] == 'E' || argv[i][0] == 'e') {
				method = METHOD_EV;
			} else {
				return -1;
			}
		} else {
			return -1;
		}
	}

	if (fiber_pathname == NULL || dst_pathname == NULL || fiber_length == 0.0f || 
		kernel_size == 0.0f || mode == -1 || method == METHOD_NONE)
		return -2;

	return 0;
}

void Inv33(double *src, double *dst)
{
	double p1 = src[0] * (src[4] * src[8] - src[5] * src[7]);
	double p2 = src[1] * (src[3] * src[8] - src[5] * src[6]);
	double p3 = src[2] * (src[3] * src[7] - src[4] * src[6]);

	double det = p1 - p2 + p3;

	double t[9];
	t[0] = src[0]; t[1] = src[3]; t[2] = src[6];
	t[3] = src[1]; t[4] = src[4]; t[5] = src[7];
	t[6] = src[2]; t[7] = src[5]; t[8] = src[8];

	dst[0] = (t[4] * t[8] - t[5] * t[7]) / det * 1.0;
	dst[1] = (t[3] * t[8] - t[5] * t[6]) / det * -1.0;
	dst[2] = (t[3] * t[7] - t[4] * t[6]) / det * 1.0;

	dst[3] = (t[1] * t[8] - t[2] * t[7]) / det * -1.0;
	dst[4] = (t[0] * t[8] - t[2] * t[6]) / det * 1.0;
	dst[5] = (t[0] * t[7] - t[1] * t[6]) / det * -1.0;

	dst[6] = (t[1] * t[5] - t[2] * t[4]) / det * 1.0;
	dst[7] = (t[0] * t[5] - t[2] * t[3]) / det * -1.0;
	dst[8] = (t[0] * t[4] - t[1] * t[3]) / det * 1.0;
}

/* CUDA functions */
int GetMaximumGflopsDeviceID(void)
{
	int gpu_num = 0;
	cutilSafeCall(cudaGetDeviceCount(&gpu_num));
	if(gpu_num > MAX_GPU_COUNT) gpu_num = MAX_GPU_COUNT;
	printf("\n\tCUDA-capable device count: %i\n", gpu_num);

	/* find maximum multiprocessors */
	int max_multiprocessors = 0;
	int max_ID = 0;
	for (int deviceID = 0; deviceID < gpu_num; ++deviceID) {
		cudaDeviceProp properties;
		cutilSafeCall(cudaGetDeviceProperties(&properties, deviceID));
		if (max_multiprocessors < properties.multiProcessorCount) {
			max_multiprocessors = properties.multiProcessorCount;
			max_ID = deviceID;
		}
	}

	printf("\tUsing No.%d GPU with %d multi processors\n\n", 
		max_ID, max_multiprocessors);

	return max_ID;
}

int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}
	
	/* read the configure file */
	ReadConfigureFile(fiber_pathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);
	
	int deviceID = GetMaximumGflopsDeviceID();
	// set the device
	cutilSafeCall(cudaSetDevice(deviceID));

	///* compute structure tensor for each voxel */
	//pTempTensor = new CDTI();
	//pTempTensor->CreateDTI(w, h, d);
	//for (int z = 0; z < d; ++z) {
	//	for (int y = 0; y < h; ++y) {
	//		for (int x = 0; x < w; ++x) {
	//			float dti[7] = {1.0f, 1.0f, 0.0f, 0.0f, 1.0f, 0.0f, 1.0f};
	//			pTempTensor->SetDTI(x, y, z, dti);
	//		}
	//	}
	//}

	//pStructureTensor = new CDTI();
	//pStructureTensor->CreateDTI(w, h, d);

	if (method == METHOD_POSITION) {
		ComputeStructureTensorMethodPosition(fiber_length, maxLength, stepSize, 
			w, h, d, scale, mode, confPathname, vertPathname, tempPathname);
	} else if (method == METHOD_FA && fa_pathname != NULL) {
		ComputeStructureTensorMethodFA(fiber_length, maxLength, stepSize, 
			w, h, d, scale, mode, confPathname, vertPathname, tempPathname, fa_pathname);
	} else if (method == METHOD_EV) {
		ComputeStructureTensorMethodEV(fiber_length, maxLength, stepSize, 
			w, h, d, scale, mode, confPathname, vertPathname, tempPathname);
	} else {
		;
	}
	
	
    /* save the result */
    //pStructureTensor->SaveDTIFile(dst_pathname);
	
	for (int z = 0; z < d; ++z) {
		SafeDeleteArray(confPathname[z]);
		SafeDeleteArray(vertPathname[z]);
	}
	SafeDeleteArray(confPathname);
	SafeDeleteArray(vertPathname);

	//SafeDelete(pStructureTensor);
	//SafeDelete(pTempTensor);
}

